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ABSTRACT 

O 
<N 

J>-j| We study the dynamical stability of planetary systems consisting of one 

hypothetical terrestrial mass planet (1 or 10M ffi ) and one massive planet 
(10 M® — 10Mj U p). We consider masses and orbits that cover the range of ob- 
served planetary system architectures (including non-zero initial eccentricities), 
determine the stability limit through N-body simulations, and compare it to the 

prl [ analytic Hill stability boundary. We show that for given masses and orbits of 

c^ | a two planet system, a single parameter, which can be calculated analytically, 

describes the Lagrange stability boundary (no ejections or exchanges) but which 
diverges significantly from the Hill stability boundary. However, we do find that 
the actual boundary is fractal, and therefore we also identify a second parameter 
which demarcates the transition from stable to unstable evolution. We show the 
portions of the habitable zones of p CrB, HD 164922, GJ 674, and HD 7924 which 
can support a terrestrial planet. These analyses clarify the stability boundaries 
in exoplanetary systems and demonstrate that, for most exoplanetary systems, 
numerical simulations of the stability of potentially habitable planets are only 
necessary over a narrow region of parameter space. Finally we also identify and 

O ' provide a catalog of known systems which can host terrestrial planets in their 



cd 



O 



X 



habitable zones. 

Subject headings: stars: planetary systems - methods: n-body simulations 

1. Introduction 

The dynamical stability of extra-solar planetary systems can constrain planet formation 
models, reveal commonalities among planetary systems and may even be used to infer the 
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existence of unseen companions. Many authors h ave studied the dynamical stability of our 



solar system and extra - solar planetary systems (see Wisdom 1982 - L Laskai 19891 ; JRasio fc Ford 
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for exam- 
ple). These investigations have revealed that planetary systems are close to dynamical insta- 
bility, illuminated the boundaries between stable and unstable configurations, and identified 
the parameter space that can support additional planets. 

From an astrobiological point of view, dynamically stable habitable zones (HZs) for 
terrestrial mass planets (0.3M® < M p < 10 M e ) are the most interesting. Classically, the 
HZ is defined as the circumstellar region in which a terrestrial m ass planet wi t h favorable 



atmospheric condit i ons can sustain liquid water on its surface (jHuangi Il959l ; iHartl 11978 
Kasting et al.lll993l ; ISelsis et al.l 120071 . but see also Barnes et al.(2009)). 



Previous work (IJones et al.ll200ll ; iMenou fc Tabachnikll2003l ; IJones et al.ll2006l ; ISandor et al 



20071 ) investigated the orbital stability of Earth- mass planets in the HZ of systems with a 
Jupiter-mass companion. In their pioneering work, IJones et al.l (120011) estimated the stability 
of four known planetary systems in the HZ of their host stars. IMenou fc Tabachniki (120031 ) 
considered the dynamical stability of 100 terrestrial mass-planets (modelled as test particles) 
in the HZs of the then-known 85 extra-solar planetary systems. From their simulations, they 
generated a tabular list of stable HZs for all observed systems. However, that study did 
not systematically consider eccentricity, is not generalizable to arbitrary planet masses, and 
relies on numerical experiments to determine stability. A similar study by I Jones et al.l (120061 ) 
also examined the stability of Earth-mass planets in the HZ. Their results indicated that 
41% of the systems in their sample had "sustained habitability" . Their simulations were 
also not generalizable and based on a large set of numerical experim ents which assumed 
the potentially habitable planet was on a circular orbit. Most recently, ISandor et al.l (120071 ) 
considered systems consisting of a giant planet with a maximum eccentricity of 0.5 and a 
terrestrial planet (modelled as a test particle initially in circular orbit) They used relative 
Lyapunov indicators and fast Lyapunov indicators to identify stable zones and generated a 
stability catalog, which can be applied to systems with mass-ratios in the range 10~ 4 — 10~ 2 
between the giant planet and the star. Although this catalog is generalizable to massive 
planets between a Saturn-mass and 10Mj up , it still assumes the terrestrial planet is on a 
circular orbit. 

These studies made great strides toward a universal definition of HZ stability. However, 
several aspects of each study could be improved, such as a systematic assessment of the 
stability of terrestrial planets on eccentric orbits, a method that eliminates the need for 
computationally expensive numerical experiments, and wide coverage of planetary masses. 



-3- 



In this investigation we address each of these points and develop a simple analytic approach 
that applies to arbitrary configurations of a giant-plus-terrestrial planetary system. 

As of March 2010, 376 extra-solar planetary systems have been detected, and the ma- 
jority (331, ~ 88%) are single planet systems. This opens up the possibility that there may 
be additional pl a nets not yet detected, in the stable regions of these systems. According 
to 



Wright et al.l (120071 ) more than 30 % of known s ingle p lanet systems show evidence for 



additional companions. Furthermore, iMarcy et al.l (j2005al ) showed tha t the distribu t ion o f 
ob served planet s rises s teeply towards smaller masses. The analyses of IWright et al.l (120071 ) 
& IMarcy et al.l ( 2005al ) suggests that many systems may have low mass planetsj. There- 
fore, maps of stable regions in known planetary systems can aid observers in their quest to 
discover more planets in known systems. 

We consider two definitions of dynamical stability: (1) Hill stability: A system is Hill- 
stable if the ordering of planets is conserved, even if the outer-most planet escapes to infinity. 
(2) Lagrange stability: In this kind of stability, every planet's motion is bounded, i.e, no 
planet escapes from the system and exchanges are forbidden. Hill st ability for a two-planet 



non-resonant sy stem can be described by an analytical expression (jMarchal fc Bozislll982 



Gladmanl 119931 ) , whereas no analytical criteria are available f or Lagrange stability so we 



inves tigate it through numerical simulations. Previous studies by lBarnes fc Greenbergl (12006 



20071 ) showed that Hill stability is a reasonable approximation to Lagrange stability in the 
case of two approximately Jupiter-mass pla nets. Part of the goa l of our pres ent work is to 
broaden the parameter space considered by iBarnes fc Greenbergj (120061 . 120071 ) . 



In this investigation, we explore the stability of hypothetical 1 M® and 10 M® planets 
in the HZ and in the presence of giant and super-Earth planets. We consider nonzero initial 
eccentricities of terrestrial planets and find that a modified version of the Hill stability 
criterion adequately describes the Lagrange stability boundary. Furthermore, we provide 
an analytical expression that identifies the Lagrange stability boundary of two-planet, non- 
resonant systems. 

Utilizing these boundaries, we provide a catalog of fractions of HZs that are Lagrange 
stable for terrestrial mass planets in all the currently known single planet systems. This 
catalog can help guide observers toward systems that can host terrestrial-size planets in 
their HZ. 



Wittenmver et al. ( 2009 ) d id a c omprehensive study of 22 planetary systems using the Hobby- Eberly 



Telescope (HET; lRamsev et al.l (J1998I )) and found no additional planets, but their study had a radial velocity 
(RV) precision of just 10 ~ 20 ms _1 , which can only detect low-mass planets in tight orbits. 
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The plan of our paper is as follows: In Section |2j we discuss the Hill and Lagrange 
stability criteria, describe our numerical methods, and present our model of the HZ. In 
Section [31 we present our results and explain how to identify the Lagrange stability boundary 
for any system with one > 10 M e planet and one < 10 M e planet. In Section HJ we apply 
our results to some of the known single planet systems. Finally, in Section [5j we summarize 
the investigation, discuss its importance for observational programs, and suggest directions 
for future research. 



2. Methodology 



According to lMarchal &: Bozisl ( 119821 ). a system is Hill stable if the following inequality 
is satisfied: 



2M 
G 2 Ml 



c 2 h > 1 + 3 4 / 3 - 



mim 2 



2/3, 



(1) 



m 3 ' (mi + ^i 2 ) 4 / 3 

where M is the total mass of the system, G is the gravitational constant, M* = mim 2 + 
m2m3 + m^rrii, c is the total angular momentum of the system, h is the total energy, mi, 
m 2 and m^ are the masses of the planets and the star, respectively. We call the left hand 
side of Eq.flT]) (3 and the right hand-side (3 crit . If (3/(3 crit > 1, then a system is definitely Hill 
stable, if not the Hill stability is unknown. 



Studies by iBarnes fc Greenberg) (120061 . 120071 ) found that for two Jupiter mass planets, 
if /5//3 cr it >1 (and no resonances are present), then the system is Lagrange stable. More- 
over, Barnes et al. ('2008a ) found that systems t end t o be packed if (3/(3 crit < 1.5 and not 



packed when (3/ /3 crit > 2. IBarnes fc Greenberg! (120071 ) pointed out that the vast majority 
of two-planet sys t ems are observed with /3//3 crit < 1.5 and hence are packed. Recently, 
Kopparapu et al.l (120091 ) proposed that the HD 47186 planetary system, with /3//3 cr i t = 6.13 
the largest value among known, non-controversial systems that have not been affected by 
tideaj, may have at least one additional (terrestrial mass) companion in the HZ between the 
two known planets. 

To determine the dynamically stable regions around single planet systems, we numer- 
ically explore the mass, semi-major axis and eccentricity space of model systems, which 
cover the range of observed of extra-solar planets. In all the models (listed in Table [1]), we 



2 See http : //xsp . astro . Washington . edu for an up to date list of /?//3 c 
multiple planet systems. 



values for the known extra-solar 
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assume that the hypothetical additional planet is either 1 M ffi or 10 M ffi and consider the 
following massive companions, (which we presume are already known to exist): (1) 10Mj up , 
(2) 5.6M jup (3) 3M jup , (4) 1.77 M jup (5) lM jup , (6) 1.86M saturn (7) 1 M saturn , (8) 56 M e (9) 
30 M e (10) 17.7 M ffi and (11) 10 M e . Most simulations assume that the host star has the 
same mass, effective temperature (T e g) and luminosity as the Sun. Orbital elements such as 
the longitude of periastron w are chosen randomly before the beginning of the simulation 
(Eq. (1) only depends weakly on them). For "known" Saturns and super-Earths, we fix 
semi-major axis a at 0.5 AU (and the HZ is exterior) or at 2 AU (the HZ interior). For 
super- Jupiter and Jupiter mass, a is fixed either at 0.25 AU or at 4 AU. These choices allow 
at least part of the HZ to be Lagrange stable. Although we choose configurations that focus 
on the HZ, the results should apply to all regions in the system. 

We explore dynamical stability by performing a large number of N-body simulations, 
each with a different initial condition. For the known planet, we keep a constant and vary 
its initial eccentricity, e, from — 0.6 in steps of 0.05. We calculate (3/(3 CT n from Eq.flT]), 
by varying the hypothetical planet's semi-major axis and initial eccentricity. In order to 
find the Lagrange stability b oundary, we perf orm numerical simulations along a particular 
/V&rit curve ; with Mercury ( iChambersI Il999l ). using the hybrid integrat or. We integrate 
each configuration for 10 7 years, long enough to identify unstable regions ( IBarnes fc Quinn 
20041 ) . The time step was small enough that energy is conserved to better than 1 part in 
10 6 . A system is considered Lagrange unstable if the semi-major axis of the terrestrial mass 
planet changes by 15% of the initial value or if the two planets come within 3.5 Hill radii of 
each otheiQ In total we ran ~ 70, 000 simulations which required ~ 35, 000 hours of CPU 
time. 



We use the definition of the "eccentric habitable zone" (EHZ; IBarnes et al.l (j2008bl )). 
which is the HZ from Selsis et al (2007), with 50% cloud cover, but assumes the orbit- 
averaged flux determines surface temperature (Williams & Pollard 2002). In other words, 
the EHZ is the range of orbits for wh ich a planet receiv es as much flux over an orbit as a 
planet on a circular orbit in the HZ of lSelsis et al.l ( 120071 ). 



3 A recent study by ICuntz fc Yeagerl ( 20091 ) notes that the Hill-radius criterion for ejection of a Earth 
mass planet around a giant planet may not be valid. Our stability maps shown here are, therefore, accurate 
to within the constraint highlighted by that study. 
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3. Results: Dynamical Stability in and around Habitable Zones 

3.1. Jupiter mass planet with hypothetical Earth mass planet 

In Figs. 1 & 2, we show representative results of our numerical simulations from the 
Jupiter mass planet with hypothetical Earth mass planet case discussed in Section [2j In 
all panels of Figs. 1 & 2, the blue squares and red triangles represent Lagrange stable 
and unstable simulations respectively, the black circle represents the "known" planet and 
the shaded green region represents the EHZ. For each case, we also plot /3//3 crit contours 
calculated from Eq. (pQ). In any given panel, as a increases, the curves change from all 
unstable (all red triangles) to all stable (all blue squares), with a transition region in between. 

We designate a particular (3/ (3 CTit contour as r s , beyond which (larger values) a hypo- 
thetical terrestrial mass planet is stable for all values of a and e, for at least 10 7 years. We 
tested r s is the first /3//3 cr it (close to the known massive planet) that is completely stable 
(only blue squares). For /3//3 crit curves below r s , all or some locations along those curves may 
be unstable; hence, r s is a conservative representation of the Lagrange stability boundary. 
Similarly, we designate r u as the largest value of /3//3 crit for which all configurations are un- 
stable. Therefore, the range t u < /5//5 crit < r s is a transition region, where the hypothetical 
planet's orbit changes from unstable (r u ) to stable (r s ). Typically this transition occurs over 
10~ 3 /3//3 cr i t . Although Figs. 1 & 2 only show curves in this transition region, we performed 
many more integrations at larger and smaller values of /3//5 crit , but exclude them from the 
plot to improve the readability. For all cases, all our simulations with /3//3 cr it > r s o,re stable, 
and all with (3/(3 crit < r u are unstable. 

These figures show that the Lagrange stability boundary significantly diverges from Hill 
stability boundary, as the eccentricity of the known Jupiter-mass planet increases. Moreover, 
r s is more or less independent (within 0.1%) of whether the Jupiter mass planet lies at 0.25 
AU or at 4 AU. If an extra-solar planetary system is known to have a Jupiter mass planet, 
then one can calculate /5//3 crit over a range of a & e, and those regions with (3/ f3 CTit > r s are 
stable. We show explicit examples of this methodology in Section 4. 

We also consider host star masses of O.3M and performed additional simulations. We 
do not show our results here, but they indicate the mass of the star does not effect stability 
boundaries. 
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3.2. Lagrange stability boundary as a function of planetary mass & 

eccentricity. 

In this section we consider the broader range of "known" planetary masses discussed in 
Section 2 and listed in Table [TJ Figures [1] & [2] show that as the eccentricity of the "known" 
planet increases, r s and r u appear to change monotonically. This trend is apparent on all our 
simulations, and suggests r s and r u may be described by an analytic function of the mass 
and eccentricity of the known planet. Therefore, instead of plotting the results from these 
models in a — e space, as shown in Figs. [U&El we identified these analytical expressions that 
relate r s and r u to mass mi and eccentricity e\ of the known massive planet. Although these 
fits were made for planets near the host star, these fits should apply in all cases, irrespective 
of it's distance from the star. In the following equations, the parameter x = log[mi], where 
??7,i is expressed in Earth masses and y — e\. The stability boundaries for systems with 
hypothetical 1 M ffi and 10 M© mass companion are: 

c 2 C4 2 y c-j , y 2 y 

7j = ci H h c 3 y + — + c 5 y + Cq - + — + c 8 y + c 9 V c w — (2) 

where j — s,u indicate stable or unstable and the coefficients for each case are given Table 
2. 

The coefficients in the above expression were obtained by finding a best fit curve to our 
model data that maximizes the R 2 statistic, 

R 2 = 1--^ (3) 



£< 



T model _ T model\'i 



where r"^ odel is the 2 ih model value of r from numerical simulations, 17 is the corresponding 
model value from the curve fit, r model is the average of all the r mo dei values and n = 572 
is the number of models (including mass, eccentricity and locations of the massive planet). 
Values close to 1 indicate a better quality of the fit. In Fig. |3l the top panels (a) & (b) 
show contour maps of r s as a function of log [mi] and e\ between model data (solid line) and 
best fit (dashed line). The R 2 values for 1 M e companion (Fig. EJ^a)) and 10M ffi companion 
(Fig. [3(b)) are 0.99 and 0.93, respectively, for r s . In both the cases, the model and the 
fit deviates when the masses of both the planets are near terrestrial mass. Therefore our 
analysis is most robust for more unequal mass planets. The residuals between the model and 
the predicted r s values are also shown in Fig. [3(c) (1M© companion) & Fig. [3(d) (10 M© 
companion). The standard deviation of these residuals is 0.0065 and 0.0257 for 1M© and 



10 M®, respectively, though the 1 M® case has an outlier which does not significantly effect 
the fit. The maximum deviation is 0.08 for 1 M® and 0.15 for 10 M® cases. 

The expression given in Eq. (J2]) can be used to identify Lagrange stable regions (/3//3crit > 
r s ) for terrestrial mass planets around stars with one known planet with e < 0.6 and may 
provide an important tool for the observers to locate these planetgj. Once a Lagrange 
stability boundary is identified, it is straightforward to calculate the range of a and e that 
is stable for a hypothetical terrestrial mass planet, using Eq.flTJ. In the next section, we 
illustrate the applicability of our method for selected observed systems. 



4. Application to observed systems 

The expressions for t s given in §3.21 can be very useful in calculating parts of HZs that 
are stable for all currently known single planet systems. In order to calculate this fraction, we 
used orbital parameters from the Exoplanet Data Explorer maintained by California Planet 
Survey consortium^, and selected all 236 single planet systems in this database with masses 
in the range 10 Mj up to 10 M© and e < 0.6. 

Table [3] lists the properties of the example systems that we consider in §4.U - §4.4] along 
with the orbital parameters of the known companions and stellar mass. The procedure to 
determine the extent of the stable region for a hypothetical 1 M® and 10 M® is as follows: (1) 
Identify the mass (mi) and eccentricity (e) of the known planet. (2) Determine r s from Eq. [2] 
with coefficients from Table 2. (3) Calculate (3//3 CT i t over the range of orbits (a and e) around 
the known planet using Eq. [TJ (4) The Lagrange stability boundary is the (3/(3 crit = t s curve. 



4.1. Rho CrB 

As an illustration of the internal Jupiter + Earth case, we consider the Rho CrB sys- 
tem. Rho CrB is a GOV star with a mass similar to the Sun, but with greater luminosity. 



Noyes et al.l ( 119971 ) discovered a Jupiter-mass planet orbiting at a distance of 0.23 AU with 
low eccentricity (e = 0.04). Since the current inner edge of the circular HZ of this star lies 
at 0.90 AU, there is a good possibility for terrestrial planets to remain stable within the 



4 Note that a more thorough exploration of mass and eccentricity parameter space may indicate regions 
of resonances on both sides of the stability. Hence, we advice caution in applying our expression in those 
regions. 



''http : //exoplanets . org/ 
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HZ. Indeed, I Jones et al.l (J200ll ) found that stable orbits may be prevalent in the present day 



circular HZ of Rho CrB for Earth mass planets. 

Fig. Hk shows the EHZ (green shaded) assuming a 50% cloud cover in the a — e space of 
Rho CrB. The Jupiter mass planet is the blue filled circle. Corresponding r s values calculated 
from Eg. (121) for 1M® companion (0.998, dashed magenta line) and 10 M e companion (1.009, 
black solid line) are also shown. These two contours represent the stable boundary beyond 
which an Earth-mass or super-Earth will remain stable for all values of a and e (cf. Figs. 
Hk anddb). The fraction of the HZ that is stable for 1 M e is 72.2% and for 10 M e is 77.0%. 
Therefore the Lagrange stable stable region is larger for a larger terrestrial planet. We 
conclude that the HZ of rho CrB can support terrestrial-mass planets, except for very high 
eccentricity (e > 0.6). 



These results are in agreement with the conclusion of lJones et al.l ( 120011 ) and lMenou Sz Tabachnik 



( 120031 ). who found that a planet with a mass equivalent to the Earth-moon system, when 
launched with e = within the HZ of Rho CrB, can remain stable for ~ 10 8 years. They 
also varied the mass of Rho CrB b up to 8.8 Mj up and still found that the HZ is stable. Our 
models also considered systems with 3Mj up , 5Mj up and 10Mj up and our results show that 
even for these high masses, if the initial eccentricity of the Earth-mass planet is less than 
0.3, then it is stable. 

To show the detectability of a 10 M® planet, we have also drawn a radial velocity (RV) 
contour of 1 ms" 1 (red curve), which indicates that a 10 M e planet in the HZ is detectable. 
A similar contour for an Earth mass planet is not shown because the precision required is 
extremely high. 



4.2. HD 164922 



Butler et al.l (120061 ) discovered a Saturn-mass planet (0.36 Mj up ) orbiting HD 164922 
with a period of 1150 days (a = 2.11 AU) and an eccentricity of 0.1. Although it has a low 
eccentricity, the uncertainty (0.14) is larger than the value itself. Therefore, the appropriate 
Saturn mass cases could legitimately use any e in the range 0.0 < e < 0.25 cases, but we use 
e = 0.1. 

Figure Hb, shows the stable regions in the EHZ (green shaded) of HD 164922, for hy- 
pothetical Earth (magenta) and super-Earth (black) planets. The Saturn-mass planet (blue 
filled circle) is also shown at 2.11 AU. About 28% of the HZ in HD 164922 is stable for a 
10 M e planet (for eccentricities <0.6), whereas for Earth mass planets only 10% of the HZ 
is stable. We again show the detection limit for a 10 M e case. 
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4.3. GJ 674 

GJ 674 is an M- dw arf star with a mass of O.35M and an effective temperature of 3600 



K. iBonfils et al.l ( 120071 ) found a 12M ffi with an orbital period and eccentricity of 4.69 days 
(a = 0.039 AU) and 0.20, respectively. Fig. Hfc shows the EHZ of GJ 674 in a — e space. Also 
shown are the known planet GJ 674 b (filled blue circle), EHZ (green shaded), and detection 
limit for an Earth- mass planet (red curve). The values of r s for 1 M e and 10 M ffi planets, from 
EqJ21 are 0.973 (magenta) and 1.0 (black), respectively. Notice that the fraction of the HZ 
that is stable for 1M© is slightly greater (79.1%) than 10M ffi planet (78.8%), which differs 
from previous systems we considered here. A similar behavior can be seen in another system 
(HD 7924) that is discussed in the next section. It seems that when the planet mass ratio 
is approaching 1, the HZ of a 10M ffi mass planet offers less stability at high eccentricities 
(> 0.6) than a 1 M e planet. But as noted in §3.2[ this analysis should be weighted with the 
fact that our fitting procedure is not as accurate for a 10 M e planet than a 1 M ffi planet. 



4.4. HD 7924 

Orbiting a K0 dwarf star at 0.057 AU, the super-Earth HD 7924 b was discovered by 



NAS A-UC Eta-Earth survey by the California Planet Search (CPS) group fjHoward et al. 



20091 ). in an effort to find planets in the mass range of 3 — 30M®. It is estimated to have 
an Msinz = 9.26M e with an eccentricity of 0.17. Fig. HJi shows r s values for hypothetical 
10M e (magenta) and lM e (black line) planets are 1.00 and 0.98, respectively. Unlike GJ 
674, where only part of the HZ is stable, around 94% of HD 7924's HZ is stable for these 
potential planets. Furthermore, we have also plotted an RV contour of 1 ms" 1 arising from 
the 10M e planet (red curve). This indicates that this planet may lie above the current 
detection threshold, and may even be in the HZ. 



Howard et al.l (120091 ) do find some additional best-fit period solutions with very high 
eccentricities (e > 0.45), but combined with a false alarm probability (FAP) of > 20%, they 
conclude that these additional signals are probably not viable planet candidates. Further 
monitoring may confirm or forbid the existence of additional planets in this system. 



4.5. Fraction of stable HZ 

For astrobiological purposes, the utility of r s is multi-fold. Not only is it useful in 
identifying stable regions within the HZ of a given system, but it can also provide (based 
on the range of a & e) what fraction of the HZ is stable. We have calculated this fraction 
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for all single planet systems in the Exoplanet Data Explorer, as of March 25 2010. The 
distribution of fractions of currently known single planet systems is shown in Fig. [5] and 
tabulated in Table |H A bimodal distribution can clearly be seen. Nearly 40% of the systems 
have more than 90% of their HZ stable and 38% of the systems have less than 10% of their 
HZ stable. Note that if we include systems with masses > 10Mj up and also e > 0.6 (which 
tend to have a ~ 1 AU (IWright et al.ll2009l )). the distribution will change and there will be 
relatively fewer HZs that are fully stable. 



5. Summary 

We have empirically determined the Lagrange stability boundary for a planetary system 
consisting of one terrestrial mass planet and one massive planet, with initial eccentricities 
less than 0.6. Our analysis shows that for two-planet systems with one terrestrial like planet 
and one more massive planet, Eq.([2]) defines Lagrange stable configurations and can be used 
to identify systems with HZs stable for terrestrial mass planets. Furthermore, in Table H] we 
provide a catalog of exoplanets, identifying the fraction of HZ that is Lagrange stable for 
terrestrial mass planets. A full version of the table is available in the electronic edition of 
the journajj. 

In order to identify stable configurations for a terrestrial planet, one can calculate a 
stability boundary (denoted as r s from Eq.([2])) for a given system (depending on the eccen- 
tricity and mass of the known planet), and calculate the range of a & e that can support 
a terrestrial planet, as shown in Section 4. For the transitional region between unstable 
to stable (r u < /?//3 crlt < r s ), a numerica l integration should be ma d e. Our results are in 
general agreement w ith previous studies (IMenou &; Tabachnikl|2003l ). (I Jones et al.ll2006l ) & 
fjSandor et all 120071 ). but crucially our approach does not (usually) require a large suite of 
N-body integrations to determine stability. 

We have only considered two-planet systems, but the possibility that the star hosts 
more, currently undetected planets is real and may change the stability boundaries outlined 
here. However, the presence of additional companions will likely decrease the size of the 
stable regions shown in this study. Therefore, those systems that have fully unstable HZs 
from our analysis will likely continue to have unstable HZs as more companions are detected 
(assuming the mass and orbital parameters of the known planet do not change with these 
additional discoveries). The discovery of an additional planet outside the HZ that destabilizes 
the HZ is also an important information. 



Updates to this catalog is available at http://gravity.psu.edu/~ravi/planets/ 
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As more extra-solar planets are discovered, the resources required to follow-up grows. 
Furthermore, as surveys push to lower planet masses, time on large telescopes is required, 
which is in limited supply. The study of exoplanets seems poised to transition to an era in 
which systems with the potential to host terrestrial mass planets in HZs will be the focus of 
surveys. With limited resources, it will be important to identify systems that can actually 
support a planet in the HZ. The parameter t s can therefore guide observers as they hunt for 
the grand prize in exoplanet research, an inhabited planet. 

Although the current work focuses on terrestrial mass planets, the same analysis can be 
applied to arbitrary configurations that cover all possible orbital parameters. Such a study 



could represent a significant improvement over the work of iBarnes fc Greenbergi (120071 ) . The 
results presented here show that /5//5 crit = 1 is not always the Lagrange stability boundary, 
as they suggested. An expansion of this research to a wider range of planetary and stellar 
masses and larger eccentricities could provide an important tool for determining the stability 
and packing of exoplanetary systems. Moreover, it could reveal an empirical relationship 
that describes the Lagrange stability boundary for two planet systems. As new planets are 
discovered in the future, the stability maps presented here will guide future research on the 
stability of extra-solar planetary systems. 
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Fig. 1.— A comparison of Hill and Lagrange stability. Colored curves (shown also in 
different linestyles) are contours of (3/(3 crit (the Hill boundary lies at f3/(3 CTit = 1). Points 
on the curve designate TV— body simulations: red points were unstable, blue stable for an 
Earth mass planet. The green shaded region represents the HZ, and the black point is the 
"known" Jupiter mass planet. The left-most curves with no stable configurations correspond 
to r u , the right-most that are fully stable represent r a . (Note that for these cases we consider 
eccentricities > 0.9 in order to identify r u ). 
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Fig. 2. — Same as Fig. dj only now the Jupiter mass planet is at 4 AU. The range of values 
of r s is approximately the same as in Fig. [TJ from 1.002 — 0.756. 
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Fig. 3. — In the top panels, we show contours of r s from numerical simulations (solid line) 
compared to the best fit (dashed line ) in log [mi] - e space, for 1 M e (panel a) and 10 M e 
(panel b) companions. The expression for the best fit is given in Eq. (J2J) with appropriate 
coefficients given in Table 1. The bottom panels show residuals between the numerical results 
and the best fit, with a standard deviation of 0.0065 for 1 M e companion (panel c) and 0.0257 
for 10 M ffi companion (panel d). 
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Fig. 4. — Comparison of Lagrange stable regions and habitable zones for four known systems. 
The magenta curves represent t s for a 1 M® planet, black 10M ffi (c.f. Fig. 3). For panels a, 
c and d, stable orbits lie to the right of these curves, but lie to the left in panel b. The red 
solid line shows the 1 ms _1 RV amplitude of a hypothetical terrestrial planet on a circular 
orbit. The green region is the HZ. 
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Fig. 5.— Distribution of fraction of stable HZ for hypothetical 1 M e planet (panel a) and 
10 M e planet (panel b), in currently observed single planet systems. Of the total systems 
(236) that we considered, nearly 40% (95) of the systems have > 90% of their HZ stable 
(peak near 1). About 38%(91) of the systems have less than 10% of their HZ stable (peak 
near 0). We do not consider systems that have planetary masses > 10Mj up or if e > 0.6. 
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Table 1. Properties of systems considered in this study. 



'Known" Planet a (AU) 
10 M jup (0.25, 4) 

5.6 M jup (0.25, 4) 



3 M jup (0.25, 4) 

1.77 M jup (0.25, 4) 

1 M jup (0.25, 4) 

1.86 M Sat (0.5, 2) 

1 M Sat (0.5, 2) 

56 M e (0.5, 2) 

30 M e (0.5, 2) 

17.7 M e (0.5, 2) 

10 M e (0.5, 2) 
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Table 2: Best fit properties for Eqf2J 




Coefficients 


lM e 


10 M e 


7"s 


r u 


r s 


Tu 


C\ 


1.0018 


1.0098 


0.9868 


1.0609 


c 2 


-0.0375 


-0.0589 


0.0024 


-0.3547 


c 3 


0.0633 


0.04196 


0.1438 


0.0105 


c 4 


0.1283 


0.1078 


0.2155 


0.6483 


c 5 


-1.0492 


-1.0139 


-1.7093 


-1.2313 


c 6 


-0.2539 


-0.1913 


-0.2485 


-0.0827 


c 7 


-0.0899 


-0.0690 


-0.1827 


-0.4456 


c 8 


-0.0316 


-0.0558 


0.1196 


-0.0279 


c 9 


0.2349 


0.1932 


1.8752 


0.9615 


ClO 


0.2067 


0.1577 


-0.0289 


0.1042 


i? 2 


0.996 


0.997 


0.931 


0.977 


a 


0.0065 


0.0061 


0.0257 


0.0141 


Max. dev. 


0.08 


0.08 


0.15 


0.05 
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Table 3. Observed parameters of example systems presented in 

System M sini a (AU) e M+(M Q ) 

Rho CrB 1.06 M jup 0.23 0.06 (± 0.028) 0.97 

HD 164922 0.36 M jup 2.11 0.05 (± 0.14) 0.94 

GJ 674 12 M e 0.039 0.20 (± 0.02) 0.35 

HD 7924 9.26 M e 0.057 0.17 (± 0.16) 0.832 
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Table 4. Lagrange stable (r s ) and unstable (r u ) boundaries, and the corresponding fraction 
of habitable zone (FHZ) stable for terrestrial mass planets in known single planet systems. 
A full version of the table is available in the electronic edition of the Astrophysical Journal. 

System mi(Mj up ) a (AU) e r s r u FHZ t s t u FHZ 

(lM ffi ) (1M 9 ) (lM e ) (10 M e ) (10 M e ) (10 M e ; 

HD 142b 1.3057 1.04292 0.26 0.9347 0.9323 0.000 0.9552 0.9320 0.000 

HD 1237 3.3748 0.49467 0.51 0.7407 0.7401 0.213 0.7549 0.7450 0.000 

HD 1461 0.0240 0.06352 0.14 0.9920 0.9780 0.976 1.0200 0.9200 0.959 

WASP-1 0.9101 0.03957 0.00 1.0022 0.9990 0.990 1.0200 0.9980 0.991 

HIP 2247 5.1232 1.33884 0.54 0.7138 0.7111 0.000 0.7490 0.7387 0.000 



